Mainly Hβ, He II, He I, [O III], [Ar IV], [N I], Si II
Copied from the 02-01-Raman-Wings-IR-Source-C notebook:
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from mpdaf.obj import Cube
import regions
import sys
from whispy import moments
from whispy import extract
sns.set_context("talk")
sns.set_color_codes()
moments.FIGPATH = Path("../../figs")
moments.SAVEPATH = Path("../../data")
We only have the coadded version of the Peter cube, so the seeing is a bit worse than in the sharp ESO cube
datapath = Path("/Users/will/Work/Muse-Hii-Data/SMC-NGC-346/")
fitsfilepath = datapath / "PeterZeidler" / "DATACUBE_FINAL_fwhm_cor.fits"
cube = Cube(str(fitsfilepath))
WARNING: UnitsWarning: '1e-20 erg/s/cm 2/Angstrom' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
Use a generic name for the cube over the broad wavelength range to reduce need for edits when we repurpose this later.
wavmin, wavmax = 4600, 5400
Only one cube available (unlike with the pipeline ESO version), so I am removing all reference to mcube
cube = cube.select_lambda(wavmin, wavmax)
# Position of mYSO
x, y = 116, 147
# position of B0 star SSN 43
xb, yb = 95, 163
# position of B0 star SSN 40
xbb, ybb = 121, 133
m = 3
boxes = {
"SSN 152+168": regions.BoundingBox(
iymin=y - m,
iymax=y + m,
ixmin=x - m,
ixmax=x + m,
),
# "BG 152+168": regions.BoundingBox(
# iymin=y-2*m, iymax=y+2*m, ixmin=x-2*m, ixmax=x+2*m,
# ),
"SSN 43": regions.BoundingBox(
iymin=yb - m,
iymax=yb + m,
ixmin=xb - m,
ixmax=xb + m,
),
# "BG 43": regions.BoundingBox(
# iymin=yb-2*m, iymax=yb+2*m, ixmin=xb-2*m, ixmax=xb+2*m,
# ),
"SSN 40": regions.BoundingBox(
iymin=ybb - m,
iymax=ybb + m,
ixmin=xbb - m,
ixmax=xbb + m,
),
# "BG 40": regions.BoundingBox(
# iymin=ybb-2*m, iymax=ybb+2*m, ixmin=xbb-2*m, ixmax=xbb+2*m,
# ),
"blue": regions.BoundingBox(iymin=75, iymax=140, ixmin=15, ixmax=40),
"red": regions.BoundingBox(iymin=200, iymax=250, ixmin=210, ixmax=300),
"magenta": regions.BoundingBox(iymin=10, iymax=50, ixmin=100, ixmax=150),
"green": regions.BoundingBox(iymin=10, iymax=100, ixmin=200, ixmax=300),
"cyan": regions.BoundingBox(iymin=170, iymax=210, ixmin=90, ixmax=120),
}
These are the ranges that the continum is fitted to. It is OK if they contain some weak lines.
wavranges = [
(4600, 4620),
(4720, 4730),
(4760, 4800),
(4937, 4947),
(4974, 4984),
(5028, 5038),
(5070, 5100),
(5120, 5145),
(5220, 5260),
(5330, 5390),
]
fig, ax = plt.subplots(figsize=(12, 12))
for label, box in boxes.items():
yslice, xslice = box.slices
spec = cube[:, yslice, xslice].mean(axis=(1, 2))
(spec / spec.data[-1]).plot(label=label)
for wavrange in wavranges:
ax.axvspan(*wavrange, alpha=0.3)
ax.legend()
ax.set(
ylim=[0.8, 2.0],
)
sns.despine()
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
There does seem a lot more structure in the "continuum" bands, as compared with the ESO cubes. Hopefully, this will not be too much of an issue.
yslice, xslice = boxes["SSN 152+168"].slices
subcube = cube[:, yslice, xslice]
contcube = extract.fit_continuum(
subcube,
wav_ranges=wavranges,
deg=6,
median=False,
)
extract.fit_continuum: row 0
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
fig, ax = plt.subplots()
subcube.mean(axis=(1, 2)).plot()
contcube.mean(axis=(1, 2)).plot()
for wavrange in wavranges:
ax.axvspan(*wavrange, alpha=0.3)
ax.set(ylim=[2000, 5000])
[(2000.0, 5000.0)]
yslice, xslice = boxes["green"].slices
subcube = cube[:, yslice, xslice]
contcube = extract.fit_continuum(
subcube,
wav_ranges=wavranges,
deg=6,
median=False,
)
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
extract.fit_continuum: row 0 extract.fit_continuum: row 10 extract.fit_continuum: row 20 extract.fit_continuum: row 30 extract.fit_continuum: row 40 extract.fit_continuum: row 50 extract.fit_continuum: row 60 extract.fit_continuum: row 70 extract.fit_continuum: row 80
fig, ax = plt.subplots()
subcube.mean(axis=(1, 2)).plot()
contcube.mean(axis=(1, 2)).plot()
for wavrange in wavranges:
ax.axvspan(*wavrange, alpha=0.3)
ax.set(ylim=[0, 180])
[(0.0, 180.0)]
So it turns out that a high order polynomial is needed because of the broad wavelength range. And it is better to use the full set of points rather than the median and that is faster too!
Calculating the median must be expensive.
contcube = extract.fit_continuum(
cube,
wav_ranges=wavranges,
deg=6,
median=False,
)
extract.fit_continuum: row 0
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/numpy/polynomial/polyutils.py:299: RuntimeWarning: divide by zero encountered in double_scalars off = (old[1]*new[0] - old[0]*new[1])/oldlen /Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/numpy/polynomial/polyutils.py:300: RuntimeWarning: divide by zero encountered in true_divide scl = newlen/oldlen /Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/numpy/polynomial/polyutils.py:368: RuntimeWarning: invalid value encountered in add return off + scl*x /Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/numpy/polynomial/chebyshev.py:1671: RankWarning: The fit may be poorly conditioned return pu._fit(chebvander, x, y, deg, rcond, full, w)
extract.fit_continuum: row 10 extract.fit_continuum: row 20 extract.fit_continuum: row 30 extract.fit_continuum: row 40 extract.fit_continuum: row 50 extract.fit_continuum: row 60 extract.fit_continuum: row 70 extract.fit_continuum: row 80 extract.fit_continuum: row 90 extract.fit_continuum: row 100 extract.fit_continuum: row 110 extract.fit_continuum: row 120 extract.fit_continuum: row 130 extract.fit_continuum: row 140 extract.fit_continuum: row 150 extract.fit_continuum: row 160 extract.fit_continuum: row 170 extract.fit_continuum: row 180 extract.fit_continuum: row 190 extract.fit_continuum: row 200 extract.fit_continuum: row 210 extract.fit_continuum: row 220 extract.fit_continuum: row 230 extract.fit_continuum: row 240 extract.fit_continuum: row 250 extract.fit_continuum: row 260 extract.fit_continuum: row 270 extract.fit_continuum: row 280 extract.fit_continuum: row 290 extract.fit_continuum: row 300 extract.fit_continuum: row 310 extract.fit_continuum: row 320
fig, ax = plt.subplots(figsize=(12, 12))
for box in boxes.values():
yslice, xslice = box.slices
spec = cube[:, yslice, xslice].mean(axis=(1, 2))
cspec = contcube[:, yslice, xslice].mean(axis=(1, 2))
spec.plot()
cspec.plot(color="k")
for wavrange in wavranges:
ax.axvspan(*wavrange, alpha=0.3)
ax.set(
ylim=[0, 5000],
)
sns.despine()
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
prefix = f"../../big-data/ngc346-PZ-{wavmin:d}-{wavmax:d}-cube"
prefix
'../../big-data/ngc346-PZ-4600-5400-cube'
csubcube = cube - contcube
cdivcube = cube / contcube
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/mpdaf/obj/arithmetic.py:124: RuntimeWarning: invalid value encountered in true_divide var /= (b_data ** 4)
csubcube.write(
f"{prefix}-contsub.fits",
savemask="nan",
)
cdivcube.write(
f"{prefix}-contdiv.fits",
savemask="nan",
)
contcube.write(
f"{prefix}-cont.fits",
savemask="nan",
)
Compare bow shock to median of cube
mboxes = {
"sw filament": regions.BoundingBox(
iymin=30,
iymax=50,
ixmin=300,
ixmax=330,
),
"bow shock": regions.BoundingBox(
iymin=165,
iymax=205,
ixmin=240,
ixmax=290,
),
"w filament": regions.BoundingBox(
iymin=100,
iymax=130,
ixmin=25,
ixmax=55,
),
"c filament": regions.BoundingBox(
iymin=195,
iymax=210,
ixmin=155,
ixmax=195,
),
}
yslice, xslice = mboxes["bow shock"].slices
bowspec_sub = csubcube[:, yslice, xslice].mean(axis=(1, 2))
avspec_sub = csubcube.mean(axis=(1, 2))
medspec_sub = csubcube.median(axis=(1, 2))
WARNING: AstropyDeprecationWarning: The slices function is deprecated and may be removed in a future version.
Use get_overlap_slices instead. [warnings]
lowspec_sub = medspec_sub.copy()
lowspec_sub.data = np.percentile(csubcube.data, 25, axis=(1, 2))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/numpy/lib/function_base.py:4650: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. arr.partition(
nobs_subcube = csubcube.data.copy()
nobs_subcube[:, yslice, xslice] = np.nan
nobs_medspec_sub = medspec_sub.copy()
nobs_medspec_sub.data = np.nanmedian(nobs_subcube, axis=(1, 2))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/numpy/lib/function_base.py:3823: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. part.partition(kth)
hi_wavs = [4641.81, 4650.00, 4685.71, 4711.37, 4740.17]
med_wavs = [4658.10, 4701.62, 4713.14, 4861.32, 4921.93, 4931.32]
fig, ax = plt.subplots(figsize=(12, 6))
wav1, wav2 = 4620, 4950
(bowspec_sub).subspec(wav1, wav2).plot(linewidth=5, alpha=1.0, color="c")
(medspec_sub).subspec(wav1, wav2).plot(color="k", linewidth=2)
(nobs_medspec_sub).subspec(wav1, wav2).plot(color="r", linewidth=2)
for wav in hi_wavs:
ax.axvline(
wav * (1.0 + 160.0 / 3e5), ymin=0.3, ymax=0.4, color="r", linewidth=3, alpha=1.0
)
for wav in med_wavs:
ax.axvline(
wav * (1.0 + 160.0 / 3e5),
ymin=0.25,
ymax=0.3,
color="b",
linewidth=3,
alpha=1.0,
)
ax.axhline(0.0, linestyle="dashed", color="r", linewidth=1)
ax.set(
ylim=[-30, 300],
)
sns.despine()
fig.savefig(moments.FIGPATH / f"ngc346-bow-shock-BZ-spec-{wav1}-{wav2}.pdf")
These are a lot harder to see than in the ESO cube, unfortunately. In following cells I will experiment with subtracting a median (or other percentile) spectrum. In the meantime, these are the original line IDs.
Left to right:
fig, ax = plt.subplots(figsize=(12, 6))
wav1, wav2 = 4620, 4950
(bowspec_sub - medspec_sub).subspec(wav1, wav2).plot(linewidth=5, alpha=1.0, color="c")
(bowspec_sub - nobs_medspec_sub).subspec(wav1, wav2).plot(
linewidth=2, color='k', alpha=1.0)
for wav in hi_wavs:
ax.axvline(
wav * (1.0 + 160.0 / 3e5), ymin=0.3, ymax=0.4, color="r", linewidth=3, alpha=1.0
)
for wav in med_wavs:
ax.axvline(
wav * (1.0 + 160.0 / 3e5),
ymin=0.25,
ymax=0.3,
color="b",
linewidth=3,
alpha=1.0,
)
ax.axhline(0.0, linestyle="dashed", color="r", linewidth=1)
ax.set(
ylim=[-10, 100],
)
sns.despine()
fig.savefig(moments.FIGPATH / f"ngc346-bow-shock-BZ-medsub-spec-{wav1}-{wav2}.pdf")
Subtracting the median gives much better results (blue line).
I had experimented with using a lower quartile than the median for subtracting, but that was worse. I also tries doing the median over only the non-bowshock parts ogf the map, but that makes hardly any difference (black line).
The long wavelength side of the strong [O III] lines is less interesting:
fig, ax = plt.subplots(figsize=(12, 6))
wav1, wav2 = 5015, 5400
(bowspec_sub - medspec_sub).subspec(wav1, wav2).plot(linewidth=5, alpha=1.0, color="c")
(bowspec_sub - nobs_medspec_sub).subspec(wav1, wav2).plot(color="k", linewidth=2)
# for wav in ariv_lines:
# ax.axvline(wav*(1.0 + 160.0/3e5), ymin=0.3, ymax=0.4, color="k", linewidth=3, alpha=0.5)
ax.axhline(0.0, linestyle="dashed", color="r", linewidth=1)
ax.set(
ylim=[-2, 20],
)
sns.despine()
He I line at 5047.74 and [Ar III] line at 5193.69
mom4686 = moments.find_moments(csubcube.select_lambda(4686, 4691))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
dlam, _, _ = cube.get_step()
cont4686 = contcube.select_lambda(4686, 4691).mean(axis=0)
ew4686 = dlam * mom4686[0] / cont4686
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/mpdaf/obj/arithmetic.py:124: RuntimeWarning: invalid value encountered in true_divide var /= (b_data ** 4)
cont4686.write(moments.SAVEPATH / "ngc346-BZ-cont-4686-mean.fits", savemask="nan")
ew4686.write(moments.SAVEPATH / "ngc346-BZ-ew-4686-mean.fits", savemask="nan")
ew4686.plot(vmin=-0.1, vmax=1, colorbar="v")
<matplotlib.image.AxesImage at 0x133601ca0>
cont4686.plot(vmin=0, vmax=3000, colorbar="v")
<matplotlib.image.AxesImage at 0x1340b53d0>
mom4686[0].mask = mom4686[0].mask | (cont4686.data > 5000)
mom4686[0].rebin(1).plot(vmin=0, vmax=300)
<matplotlib.image.AxesImage at 0x13417bf40>
So we can still see the bow shock, but the image is not as clear as with the ESO cube. The horizontal banding artifacts are worse.
mom_pars_4686 = dict(
restwav=4685.68,
irange=[-10, 1000],
vrange=[-250, 550],
srange=[0, 600],
)
plot_pars_4686 = dict(
ilabel="He II",
label="4686",
flabel="ngc346-heii",
**mom_pars_4686,
)
g = moments.moments_corner_plot(
mom4686,
rebin=1,
**plot_pars_4686,
)
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:119: RuntimeWarning: invalid value encountered in log10
f"log10 I({label})": np.log10(mom0.data[~m]),
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:127: RuntimeWarning: invalid value encountered in log10
f"log10 I({label})": np.log10(irange),
g = moments.moments_corner_plot(
mom4686,
rebin=4,
**plot_pars_4686,
)
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:119: RuntimeWarning: invalid value encountered in log10
f"log10 I({label})": np.log10(mom0.data[~m]),
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:127: RuntimeWarning: invalid value encountered in log10
f"log10 I({label})": np.log10(irange),
Unfortunately, the velocity spread does not decrease with binning, although the sigma spread does.
(3e5 * (mom4686[1] / 4685.68 - 1.0)).rebin(1).plot(
vmin=100,
vmax=220,
cmap="seismic",
colorbar="v",
)
<matplotlib.image.AxesImage at 0x1353f49d0>
This velocity map has serious problems with the horizontal banding.
map4686 = mom4686[0].copy()
map4686.mask = map4686.mask | (map4686.data < -100)
map4686.rebin(1).plot(
vmin=0,
vmax=500,
)
<matplotlib.image.AxesImage at 0x135afbac0>
map4686.rebin(4).plot(
vmin=0,
vmax=500,
)
<matplotlib.image.AxesImage at 0x134accc10>
mom_pars_4686 = dict(
restwav=4685.68,
irange=None,
vrange=None,
srange=None,
)
moments.save_moments_to_fits(
mom4686,
label="4686",
flabel="ngc346-PZ-heii",
**mom_pars_4686,
)
# ## Extract [Ar IV] line
mom4740 = moments.find_moments(csubcube.select_lambda(4740, 4746))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
mom4740[0].plot(vmin=100.0, vmax=400.0)
<matplotlib.image.AxesImage at 0x134a66f70>
In this case, the image is as good or better than the ESO one. Note that the minimum value is 100 since there seems to be diffuse emission over the entire map. In the ESO case, this was eliminated by definition.
mom_pars_4740 = dict(
restwav=4740.17,
irange=[50, 1000],
vrange=[100, 300],
srange=[50, 150],
)
plot_pars_4740 = dict(
ilabel="[Ar IV]",
label="4740",
flabel="ngc346-PZ-ariv",
**mom_pars_4740,
)
g = moments.moments_corner_plot(
mom4740,
rebin=4,
**plot_pars_4740,
)
These moment plots look pretty good. The velocity is well-contrained to 170 km/s in the brightest parts. It looks to be slightly redder in the fainter parts (175 to 180 km/s).
moments.save_moments_to_fits(
mom4740,
label="4740",
flabel="ngc346-PZ-ariv",
restwav=4740.17,
)
(3e5 * (mom4740[1] / 4740.17 - 1.0)).rebin(8).plot(
vmin=150.0, vmax=190.0,
cmap="seismic",
colorbar="v",
)
<matplotlib.image.AxesImage at 0x133a868b0>
This is the velocity map, which confirms that the BS seems slightly blue-shifted compared with the more diffuse [Ar IV] emission.
And the other [Ar IV] line. This is contaminated by He I, which we will have to sort out at some point.
mom4711 = moments.find_moments(csubcube.select_lambda(4711, 4717))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
mom4711[0].plot(vmin=100.0, vmax=500.0)
<matplotlib.image.AxesImage at 0x134c5ad00>
(mom4711[0].rebin(2) / mom4740[0].rebin(2)).plot(
vmin=0.5,
vmax=2.5,
colorbar="v",
cmap="magma",
)
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/mpdaf/obj/arithmetic.py:124: RuntimeWarning: invalid value encountered in true_divide var /= (b_data ** 4)
<matplotlib.image.AxesImage at 0x135a74fd0>
So despite the contamination, the ratio is roughly constant at around 1.5 in the bright parts.
moments.save_moments_to_fits(
mom4711,
label="4711",
flabel="ngc346-PZ-ariv",
restwav=4711.37,
)
So far, we have been looking at median-subtracted spectra, but have not applied the same technique to the whole cube. So I am going to try that, to see if it improves the moment maps.
cmsubcube = csubcube.copy()
cmsubcube.data -= medspec_sub.data[:, None, None]
mom4740m = moments.find_moments(cmsubcube.select_lambda(4740, 4746))
mom4686m = moments.find_moments(cmsubcube.select_lambda(4686, 4691))
mom4711m = moments.find_moments(cmsubcube.select_lambda(4711, 4717))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data) /Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data) /Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
fig, axes = plt.subplots(3, 2, figsize=(8, 15))
mom4740[0].plot(vmin=100, vmax=400.0, ax=axes[0, 0])
mom4740m[0].plot(vmin=-50.0, vmax=250.0, ax=axes[0, 1])
mom4711[0].plot(vmin=100, vmax=500.0, ax=axes[1, 0])
mom4711m[0].plot(vmin=-50.0, vmax=350.0, ax=axes[1, 1])
mom4686[0].plot(vmin=0, vmax=200.0, ax=axes[2, 0])
mom4686m[0].plot(vmin=0.0, vmax=200.0, ax=axes[2, 1])
<matplotlib.image.AxesImage at 0x135917610>
So the answer is that it makes pretty much zero difference, apart from shifting the intensity baseline in some cases.
Verdict: Not worth pursuing for these lines, but see below where we find it makes a big improvement for [Ar III].
fig, ax = plt.subplots(figsize=(12, 6))
wav1, wav2 = 5170, 5220
(bowspec_sub - medspec_sub).subspec(wav1, wav2).plot(linewidth=5, alpha=1.0, color="c")
(bowspec_sub - nobs_medspec_sub).subspec(wav1, wav2).plot(color="k", linewidth=2)
# for wav in ariv_lines:
# ax.axvline(wav*(1.0 + 160.0/3e5), ymin=0.3, ymax=0.4, color="k", linewidth=3, alpha=0.5)
ax.axhline(0.0, linestyle="dashed", color="r", linewidth=1)
ax.set(
ylim=[-2, 5],
)
sns.despine()
mom5192 = moments.find_moments(cmsubcube.select_lambda(5191, 5198))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
mom5192[0].rebin(2).plot(vmin=-20, vmax=50)
<matplotlib.image.AxesImage at 0x136b7cf40>
When I used the normal cube, there is nothing there at all. BUT, with the median-subtracted cube I get something very similar to what I got from the ESO pipeline cube.
This has the classic distribution seen in medium ionization lines, where the brightest part is in the lower left.
We do have bad artefacts in vertical strips on the left and right edges though.
moments.save_moments_to_fits(
mom5192,
label="5192",
flabel="ngc346-PZ-ariii",
restwav=5191.82,
)
fig, ax = plt.subplots(figsize=(12, 6))
wav1, wav2 = 4900, 4940
(bowspec_sub - medspec_sub).subspec(wav1, wav2).plot(linewidth=5, alpha=1.0, color="c")
(medspec_sub).subspec(wav1, wav2).plot(color="k", linewidth=2)
# for wav in ariv_lines:
# ax.axvline(wav*(1.0 + 160.0/3e5), ymin=0.3, ymax=0.4, color="k", linewidth=3, alpha=0.5)
ax.axhline(0.0, linestyle="dashed", color="r", linewidth=1)
ax.set(
ylim=[-50, 130],
)
sns.despine()
mom4922 = moments.find_moments(csubcube.select_lambda(4921, 4929))
mom4922m = moments.find_moments(cmsubcube.select_lambda(4921, 4929))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data) /Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
mom4922[0].rebin(1).plot(vmin=0.0, vmax=400.0)
<matplotlib.image.AxesImage at 0x135686400>
mom4922m[0].rebin(1).plot(vmin=-50.0, vmax=350.0)
<matplotlib.image.AxesImage at 0x1377e1b80>
This one is good enough signal to see the increase in brightness at the bow shock.
moments.save_moments_to_fits(
mom4922,
label="4922",
flabel="ngc346-PZ-hei",
restwav=4921.93,
)
fig, ax = plt.subplots(figsize=(12, 6))
wav1, wav2 = 5010, 5070
(bowspec_sub - medspec_sub).subspec(wav1, wav2).plot(linewidth=5, alpha=1.0, color="c")
(medspec_sub).subspec(wav1, wav2).plot(color="k", linewidth=2)
# for wav in ariv_lines:
# ax.axvline(wav*(1.0 + 160.0/3e5), ymin=0.3, ymax=0.4, color="k", linewidth=3, alpha=0.5)
ax.axhline(0.0, linestyle="dashed", color="r", linewidth=1)
ax.set(
ylim=[-60, 100],
)
sns.despine()
mom5048 = moments.find_moments(csubcube.select_lambda(5047, 5055))
mom5048m = moments.find_moments(cmsubcube.select_lambda(5047, 5055))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data) /Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
mom5048[0].rebin(2).plot(vmin=70.0, vmax=140.0)
<matplotlib.image.AxesImage at 0x136c403d0>
mom5048m[0].rebin(2).plot(vmin=0.0, vmax=70.0)
<matplotlib.image.AxesImage at 0x137ff4d60>
moments.save_moments_to_fits(
mom5048,
label="5048",
flabel="ngc346-PZ-hei",
restwav=5047.74,
)
mom4861 = moments.find_moments(csubcube.select_lambda(4859, 4869))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
Overview before correcting the sky:
fig, axes = plt.subplots(
2,
2,
figsize=(10, 10),
sharex=True,
sharey=True,
)
imap = mom4861[0].copy()
vmap = 3e5 * (mom4861[1] / 4861.32 - 1.0)
smap = 3e5 * (mom4861[2] / 4861.32)
# m = imap.data > 10.
# vmap.mask = vmap.mask | (~m)
# smap.mask = smap.mask | (~m)
imap.rebin(1).plot(
vmin=1e4,
vmax=3e5,
cmap="turbo",
ax=axes[0, 0],
)
vmap.rebin(1).plot(
vmin=140,
vmax=190,
cmap="seismic",
ax=axes[0, 1],
)
smap.rebin(1).plot(
vmin=40,
vmax=100,
cmap="magma",
ax=axes[1, 0],
)
imap.rebin(1).plot(
vmin=-5000,
vmax=0,
cmap="viridis",
ax=axes[1, 1],
)
bg_4861 = contcube.select_lambda(4859, 4869).mean(axis=0)
axes[1, 1].contour(
bg_4861.data,
levels=[1000],
colors="r",
)
fig.tight_layout()
Make a 30 Å window:
core_4861 = csubcube.select_lambda(4849, 4879)
mom4861 = moments.find_moments(core_4861.select_lambda(4859, 4869))
/Users/will/miniconda3/envs/py39/lib/python3.9/site-packages/whispy/moments.py:43: RuntimeWarning: invalid value encountered in sqrt mom2.data = np.sqrt(mom2.data)
That is looking great now. The first moment is very similat to Hα, which is encouraging.
mom_pars = dict(
restwav=4861.32,
irange=[8.0e3, 3e5],
vrange=[150, 180],
srange=[30, 100],
)
moments.save_moments_to_fits(
mom4861,
label="4861",
flabel="ngc346-PZ-hi",
**mom_pars,
)
plot_pars = dict(
ilabel="H I",
label="4861",
flabel="ngc346-hi",
**mom_pars,
)
g = moments.moments_corner_plot(
mom4861,
rebin=1,
**plot_pars,
)
This is extremely similar to the results for Ha, which is encouraging. We see the same 3 clumps on the I-V plane, and with very similar velocities.
The sigma is a bit larger than Ha, but that is to be expected since the instrumental broadening is larger.